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The HeH"*" cation undergoes dissociative recombination witli a free electron to produce neutral He 
and H fragments. We present calculations using ab initio quantum defects and Fane's rovibrational 
frame transformation technique, along with the methodology of Ref. to obtain the recombina- 
tion rate both in the low-energy (1-300 meV) and high-energy (ca. 0.6 hartree) regions. We obtain 
very good agreement with experimental results, demonstrating that this relatively simple method 
is able to reproduce observed rates for both indirect dissociative recombination, driven by rovibra- 
tionally autoionizing states in the low-energy region, and direct dissociative recombination, driven 
by electronically autoionizing Rydberg states attached to higher-energy excited cation channels. 
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I. INTRODUCTION 

The dissociative recombination (DR) 0, [1, 0] reaction 
HeH+ H- c~ ^ He + H has received considerable theoreti- 
cal and experimental interest. HeH+ ion chemistry is im- 
portant to the understandingof the composition of inter- 
stellar space. Observations [5| have failed to show much 
HeH+ present in the interstellar medium, and were in 
contradiction to the early prediction that this species 
would be abundant. These observations indicated that 
the low-energy DR cross section of HeH+ + e~ is indeed 
significant, though there is no valence electronic state of 
the neutral that provides a mechanism for the high rate 
at low collision energy. 

Thus, much theoretical and experimental interest has 
been focused on this process in the intervening years. Ex- 
periments 0, i, S [ia, [H El, [3 show peaks in the 
DR cross section in both the low-energy region (ca. 1- 
300 meV) and the high-energy region around 20 eV. It is 
the low-energy peak that is responsible for the destruc- 
tion of HeH+ in interstellar space. 

The low and high energy peaks are caused by two dif- 
ferent mechanisms, the indirect and direct, respectively. 
The indirect mechanism has proved more difficult to treat 
theoretically, and the initial underestimate Q of the low- 
energy DR rate was based upon a consideration of the 
direct mechanism only. 

The presence of the two mechanisms, and the impor- 
tance of this process to the understanding of interstellar 
chemistry, makes the DR of HeH"'" a prime target for 
theory. Several treatments 0, [ll, [H [H, [13, [H, [H 
have been presented, and these have had varying suc- 
cess at reproducing experimental DR rates. The study 
of Orel et al. [15[ reproduced the magnitude and size of 
the high-energy peak, but did not accurately reproduce 
the shoulder on the high-energy side. Takagi [H, [l^ has 
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FIG. 1: Schematic of HeH+ DR. The ground rovibrational 
state of the cation is plotted as dots; the cation curves, as solid 
lines; and Rydberg curves supported by the lowest and high- 
est cation curves, as dotted lines. Typical incident electron 
energies for the low-energy and high-energy DR are labeled 
Eo and Ei, respectively. 



obtained good agreement for the low-energy peak, but 
did not reproduce all the peaks and dips in the cross 
section perfectly. 

As schematic of the process is shown in Fig. [TJ Indi- 
rect dr[1] in the low-energy region is driven by transi- 
tions to rovibrationally autoionizing Rydberg states sup- 
ported by the lowest cation curve. Such states are de- 
scribed as a Rydberg electron attached to a rotationally 
or vibrationally excited, but electronically ground-state, 
IIeII+ cation core. The DR process is driven by nona- 
diabatic coupling to and among these Rydberg states. 
The DR cross section exhibits peaks corresponding to 
the locations of the autoionizing resonances, which serve 
as doorways. 

The direct DR process [2^ [2l| is driven by a transfer 
of energy among the electronic degrees of freedom only: 
an incident electron with sufficient energy may excite the 
ground-state (IS) HeH^ cation core, such that the in- 
cident electron becomes trapped in a Rydberg state sup- 
ported by a dissociative a cation core. This is a vertical 
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transition, labeled by an arrow in Fig. [1] and nonadia- 
batic coupling is unimportant to the total DR rate. Once 
the electronic transition to the metastable Rydberg state 
has occurred, the molecule follows the dissociative Ryd- 
berg potential energy surface toward breakup. Our stud- 
ies indicate that most of the high-energy DR to HeH+ 
is accounted for by the Rydberg scries converging to the 
upper singlet curve. 



II. THEORY 

Several theoretical methodologies have found rea- 
sonable success in treating the dissociative recombina- 
tion of HeH+, but there is still room for improve- 
ment. Treatments[l^ [l^, using the formalism of 
O'Mallevpll. in which the DR rate is obtained by 
calculating potential energy curves of the autoioniz- 
ing Rydberg states of HeH, may be contrasted with 
treatments [isl . [iqI u sin g rn ultichannel quantum de- 
fect theory (MQDT)|22, [H,]!!, [H, [H in which the only 
potential energy curves explicitly included are those of 
the cation, or those in which both cation and neutral 
curves are used 0, [131 . Scattering calculations using vi- 
brational close coupling were used in Ref. [2^ . 

The O'Malley treatment was used by Orel et aLp^ 
to calculate the high energy peak; those authors repro- 
duced the magnitude and basic shape of the high-energy 
peak obtained from the 1993 experiment of Sundstrom 
et al. [§] by including a total of eight doubly excited Ry- 
dberg states (six S and two 11). To calculate the low- 
energy peak under the O'Malley framework, Larson et 
aL[16l| required nonadiabatic coupling matrix elements 
between the ground and E-core Rydberg states of HcH+ ; 
their study underestimated the magnitude of the exper- 
imental result by a factor of ten to 100. 

Treatments using multichannel quantum defect the- 
ory can be considerably simpler than those using the 
O'Malley treatment because they do not necessarily re- 
quire the enumeration of the Rydberg states contributing 
to the process. Such treatments include ours and those of 
Takagi [l^, [l^ and Gubcrman . The main difference 
between our method and those of Takagi and Gubmer- 
man is that ours includes the DR process in the closed 
channel space of the MQDT calculation, whereas the oth- 
ers include DR channels in the open channel space. 

Guberman's pioneering 1994 paper employed the 
hybrid technique [131 incorporating both dissociative 
neutral and cation curves. In this method one typically 
has dissociative valence states and MQDT S-matrices 
corresponding to the Ryberg series, and one proceeds 
by coupling these in the first Born approximation. The 
coupling can be purely electronic or can arise from nona- 
diabatic coupling. 

Nonadiabatic coupling was employed in Guberman's 
treatment. He defined MQDT S-matrices using the curve 
of the D state of neutral HeH. He treated the C state, 
which lies below the D state, as an open dissociative chan- 



nel. The C and D states avoid and the nonadiabatic cou- 
pling between them is included as the discrete-continuum 
coupling as in Ref. p7j . 

Guberman's technique has several advantages, notably 
that it incorporates the correct asymptotic form of the 
wavefunction. However, it requires the identification of 
the dominant final electronic state channel of the DR 
process before doing the calculation, and thus may be 
difficult to apply to a general problem. 

Our method of calculating PR rates using MQDT is 
similar to Takagi's method [isl. [Tgj in that it incorpo- 
rates a vibrational frame transformation [H, includ- 
ing both bound and continuum vibrational wavefunc- 
tions of the HeH+ cation. Our methods differ in the 
choice and treatment of the continuum vibrational states. 
Takagi employs standing-wave boundary conditions (box 
states), and chooses a subset of these states to repre- 
sent the DR channels, interspersed among others that are 
chosen to be closed. We employ outgoing-wave bound- 
ary conditions[l|, and consider such continuum channels 
open or closed according to the real part of their complex- 
valued energies. Thus, in our method, the DR channels 
are not explicitly included in the calculated s-matrix, 
whereas in Takagi and Gubmerman's treatments, they 
are. 

Whereas Guberman's calculation incorporates the cor- 
rect asymptotic form of the wavefunction, but is dif- 
ficult to implement for a general system, our method 
and that of Takagi incorporate unphysical boundary con- 
ditions and are more easy to implement. In Takagi's 
treatment, some box continuum states are chosen to be 
open and some closed, and these are interspersed among 
each other. Therefore, in some box states the bound- 
state boundary condition is imposed for the electron, and 
for others, they are not. The calculation has Rydberg 
states attached to closed box states and this situation is 
clearly unphysical. In our calculation, we have wavefunc- 
tions with complex energy, asymptotically increasing at 
large bond length, which is also clearly unphysical. How- 
ever, in the limit that the discrctized continuum states 
have a negligibly small imaginary component to their en- 
ergy, our wavefunction corresponds to the correct physi- 
cal scattering wavefunction having outgoing wave bound- 
ary conditions in the diatomic coordinate. 

At present, our method has not been applied to the 
problem of the final-state distribution of the electronic 
states of the fragment atoms. Applying the methods of 
Guberman or Takagi to this problem is more straightfor- 
ward than it is for our method. 



A. Outgoing wave basis functions 

The original outline of our methodfH as well as sub- 
sequent calculations of DR rates for physical systems [l^ 
HOjUH used the technique of Siegert pseudostatesfs^ to 
define the outgoing-wave cation vibrational basis. How- 
ever, there are a few difficulties with the formal theory 
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presented in these papers that we would like to address. 
The Siegert-state basis is orthonormal with respect to 
integration plus a surface term. 



dR Xc.{R)xp{R) +1 



■ Xa{Ro)Xl3{Ro) 



ka + kp 



= Saf: 
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where k^ is the Siegert pseudostate wavenumber eigen- 
value for state number a. By analogy, the frame- 
transformed S-matrix has been written (see, for example, 
Eq.(6) in Rcf. [ij) 
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However, questions may be raised as to whether this ad 
hoc expression represents the proper transformation of 
the S-matrix from the R basis to the x basis. For in- 
stance, this equation does not preserve the eigenphases 
of the MQDT s-matrix, and is not invertible. 

The orthonormality properties of the Siegert pseu- 
dostates are well known j32|. and it is not hard to de- 
rive a more appropriate expression for the transformed 
S-matrix using Siegert pseudostates. However, in the 
present paper we choose to use a different method to 
enforce outgoing-wave boundary conditions in the vibra- 
tional basis. 

Complex absorbing potentials (CAP) |33l. [34 ] . 13511 or 
exterior complex scahng (ECS) [11, [13, EE El SI 

are alternative methods of defining outgoing-wave states. 
The set of eigenvectors obtained from these methods obey 
a simpler orthonormality relationship than do the Siegert 
pseudostates: they may be chosen orthonormal with re- 
spect to the inner product in which neither the bra nor 
ket is complex-conjugated, called the C-norm, 
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Unlike the orthonormality relationship of the Siegert 
states, the C-norm is an inner product and defines the 
ECS or CAP states as basis vectors of a vector space. 
The frame transformation of the fixed-nuclei S-matrix to 
the ECS or CAP basis uses this inner product and is 
straightforward: 



(4) 
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This transformation preserves the eigenphases exactly in 
the complete-basis limit. We then apply the rotational 
frame transformation of Ref. [l^l . 

Previous applications of the method outlined in Ref. [l| 
have been on indirect DR in which electron-impact dis- 
sociation is not energetically open. Thus, the outgoing- 
wave vibrational states have corresponded to closed chan- 
nels in these studies. In contrast, for the present analysis 
of the high-energy peak in HeH+ DR, electron-impact 



dissociation to H+ -I- He is energetically open (that to H 
-|- He^ niay be) and we have outgoing-wave vibrational 
functions as open channels. 

It is worth mentioning that in cither case, the asymp- 
totic form of the wavefunction is unphysical, because of 
the presence of outgoing wave basis functions with com- 
plex energy that are exponentially increasing in mag- 
nitude. For the frame-transformed S-matrix in the c- 
normed vibrational basis, each column of which corre- 
sponds to an outgoing- wave scattering state vI>+(*), 

= Xa(i?)VJfc«0 +E'^"/3X/3(^)/^?',(fc30 , (5) 

the proper statement of unitarity easily is found not to 
be 



but 



where 
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While it is easy to prove that the frame-transformed and 
channel-closed S-matrix with bound states is unitary in 
the complete basis limit, we have no proof that our con- 
structed S-matrix is guaranteed to be subunitary, though 
we have found this to be the case in the applications of 
this theory to date. 

Thus, the outline of the calculation is as follows. We 
calculate fixed-nuclei S-matrices sim,vmi{R) and trans- 
form them to the ECS basis via Eq.(|l]). We then ap- 
ply the rotational frame transformation of Chang and 
Fano [131 . We finally apply the channel-closing formula. 
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where the subscript c and o denote the closed and open 
channel subblocks of the MQDT S-matrix S", a is now a 
collective index including angular momentum quantum 
numbers, and wc introduce the notation 5^ for the phys- 
ical S-matrix. The DR cross section is then obtained from 
the unitary defect of the frame-transformed S-matrix, 



(10) 



,37 



We Boltzmann average and convolute with respect to the 
parallel and transverse spreads in the incident electron 
energy in the experiment, as described in Rcf. [3l|. 
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FIG. 2; Quantum eigendefects calculated at zero incident energy in the ground state cation channel. 



III. CALCULATION OF AB INITIO QUANTUM 
DEFECTS AND CATION ROVIBRATIONAL 
STATES 

We calculate quantum defects using the polyatomic 
UK R-matrix program [i^, which is based on the 
Swedish-Molecule and Alchemy suites of quantum chem- 
istry codes [13]. We include s,p, d, and /-wave basis func- 
tions (only s, p, and d for the high-energy peak) for the 
scattered electron in an R-matrix box of radius 14 bohr, 
calculated with the program GTOBASjH]. We use the 
cc-pvtz basis set. For calculations on the low-energy re- 
gion, we include the lone la SCF orbital of HcH+ and an 
additional six a and two pairs of tt virtual orbitals for the 
target space. The three lowest cation target states, 1 and 
2 and '^S, are included in the scattering calculation 
and defined by full CI in the target space. Penetration 
terms defined by full CI of three electrons in the target 
space are included in the scattering wavcfunction. 

For the high-energy DR calculations, the HeH+ target 
states are represented by full CI in the orbital space of a 
CASSCF calculation on the ground state, performed with 
the molpro quantum chemistry package |46|. Following 
Orel et aZ.[l5|, we use a total of four a orbitals and one 
TT orbital. We include the first three cation states and 
use an R-matrix radius of eight bohr. Penetration terms 
with full CI of three electrons are again included in the 
scattering wavefunction. 

Quantum defects at zero incident energy in the ground- 
state cation channel are plotted in Fig. [51 

To perform the frame-transformation calculation to 
obtain the low-energy DR rate, we interpolate the calcu- 
lated fixed-nuclei s-matrices over a range of i?=0.75 to 4.5 
bohr, and calculate outgoing-wave ECS vibrational states 
using the finite-element Gauss-Lobatto discrete variable 
representation [4^ . We used 16-point quadrature and 
seven elements, scaling the R coordinate at an angle of 
^TT on the sixth element (the calculation is almost fully 
converged at 10th order quadrature). We included 42 
vibrational wavefunctions; the calculations of Ref. [48| 
found 12 bound states, and our vibrational states in- 
clude seven with real energies and ten with imaginary 
energy component less than 0.0001 hartree in magnitude. 
We use the ground-state cation potential energy surface 
of Coxon and IIajigcorgiou[49| and use modified atomic 
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FIG. 3: (Color online) Unconvolved cross section calculated 
for the low-energy DR of the ground rovibrational state of 
^HeH. 
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where m i , is the mass of a helium nucleus plus one 
and a half electrons and m i , is the mass of a proton 

H2+ ^ 

plus a half electron. The results are slightly sensitive to 
such fine-tuning (which approximately accounts for nona- 
diabatic coupling) and we find this choice to produce the 
best agreement with prior calculations of HeII+ rovibra- 
tional spectra. For transitions up to = 3 or j = 3 
we find good agreement with the results of Bishop and 
Cheung[50|; with errors less than 21 cm~^. 



IV. RESULTS 

Results for the low-energy calculation on IIcII+ are 
plotted in Figs. [3l HI and[5l Figure [3] shows the uncon- 
volved cross section for DR of the ground rovibrational 
state of ^HeH. To compare with experiment, we also con- 
volve our results using a transverse spread of lOmeV 
and a parallel spread of O.lmeV in the incident electron 
ener gy, to account for the experimental parameters of 
Ref. [ii]. 
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FIG. 4: (Color online) Cross section calculated in the vicinity 
of the j = 1 threshold at 8.3meV. Unconvolved physical DR 
cross section, red line; physical DR cross section convolved 
with a gaussian, purple line; unphysical cross section calcu- 
lated by artificially opening the j = 1 channel, black line. 
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Inelastic scattering leads to rotational or vibrational 
excitation of the cation and is allowed starting at an 
energy of 8.3meV, where the first rotational threshold 
(j = I) lies. As shown in Fig. [31 at that threshold there is 
a clear step down in the DR cross section, and the inelas- 
tic cross section assumes a value near the averaged DR 
rate below threshold. This finding confirms the basic pic- 
ture of the indirect mechanism, and parallels our findings 
on the DR of LiH+ [33,[3l|. Below the j = 1 threshold, 
the only inelastic process open is dissociative recombi- 
nation; above it, rotational excitation of the cation is 
allowed, and this process competes with DR and in fact 
dominates it, decreasing the DR cross section. 

Thus, for energies below the j = 1 threshold of 8.3meV, 
the DR is seen to be driven by capture into Rydberg 
states supported by the j ~ 1 rotational state of the 
cation. These resonances are visible in Fig. [3] and have 
a sharp cut-off at the 8.3meV threshold. Another forest 
of Rydberg states terminates at the second excited rota- 
tional threshold, j = 2, at 24.9meV; at this energy, there 
is a dip visible in the convolved DR cross section, but 
a series of broad resonances then appears and increases 
the cross section. The inelastic cross section has only 
a barely visible step up at the second excited rotational 
threshold, showing that rotational excitation to j = 1 is 
much stronger than that to j = 2. 

A close-up view of the calculated cross section for DR 
of the ground rovibrational state of '^Hell is shown in 
Fig. m In this figure we also plot the result of calcula- 
tions in which the first excited rotational channel {j = 1) 
is left open even below its threshold. Such a calculation 
elminates the Rydberg series converging to the J = 1 
threshold and reveals the presence of resonances attached 
to higher channel thresholds. In this figure we plot the 
unconvolved physical DR cross section; the physical DR 
cross section convolved with a gaussian to show its aver- 
age value; and the unphysical, artificially-opened j = 1 
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FIG. 5: (Color online) Low-energy DR rate coefficient com- 
pared with the experiment of Tanabe et al.^^ and the the- 
oretical results of Takagi [ill. Guberman[§j, and Sarpal et 



channel result. Note that, as in Refs.[30|, ISlJ, there is an 
enhancement of the DR rate at a complex resonance, i.e., 
a resonance involving a lower principal quantum num- 
ber state embedded in the very high Rydberg series con- 
verging to an excited cation state threshold. In contrast 
to Rcfs.fsol. [3]| . however, the complex resonance mech- 
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anism is here seen to be supported by rotationally au- 
toionizing resonances, rather than vibrationally autoion- 
izing ones. This suggests that some of the ideas shown in 
Refs.fsol. [sij . concerning the indirect mechanism of DR 
being dominated by complex resonances, may also be 
relevant for rotationally-autoionizing states in addition 
to the originally-postulated situation of vibrationally- 
autoionizing Rydbcrg states. It is beyond the scope of 
this study to examine this point in greater detail, but 
this point is worth exploring in future research. 

The results presented above analyze the DR of rota- 
tionally and vibrationally cold, ground-state ''HeH. To 
compare with current experiments, we must account 
for the nonzero temperature of the ion source. We 
Boltzmann-average over the population of rovibrational 
states at 800° K, to account for the experimental pa- 



rameters of Ref. [IJ] , and show the results for the various 



isotopologucs in Fig. 
results of Gubcrman 



We compare with the theoretical 
, Sarpal et a/. HI], and Takagi[li 
and it seems that the present method has produced the 
results most closely matching the experimental data of 
Tanabe[ll|. 

For calculating the high-energy DR peak, several de- 
grees of complication may be included in the calculation, 
but we find that a very simple treatment is sufhcient to 
approximate the experimental results. We include nei- 
ther the energy dependence nor the R-dependence of the 
fixed-nuclei S-matrix, evaluating it at 0.6 hartree and 
l.Bap, roughly at the center of the experimentally ob- 
served peak and ground-state vibrational wavefunction. 

Results for the high-energy DR calculation arc shown 
in Fig. [51 We compare our MQDT results to the prior 
theoretical results of Ref. and the experiments of 
Refs. and We also calculate a result using the 

O'Malley framework [2l[ in which we first find the reso- 
nance positions and widths as a function of nuclear geom- 
etry. We do so by applying the MQDT channel-closing 
formula, Eq.®, to the calculated MQDT S-matriccs to 
obtain the "physical" S-matrix S^{E) in the complex- 



energy plane. We then locate its poles E^^ — i^-^-- We 
find six resonances, i = 1 — 6, all of which have sigma 
symmetry. Using the approximate formula [2l| valid for 
small widths, 



r«(i?.(i?)) 



E 



E 



(12) 



where the widths F*-*', resonance energies -E'j^'', and the 
ground-state vibrational wavefunction ^' are evaluated 
at a geometry Ri {E) consistent with a vertical transition 
such that Eo+E = E^^{Ri{E)), with Eq the initial state 
energy, we obtain a result consistent with our MQDT re- 
sult. In the lower panel of Fig. [6] we show the contribu- 
tions for each of the six sigma resonances to the total. 
Orel et aL[l5| found two 11 resonances having significant 
contribution to the total cross section, and their six S 
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FIG. 6: Top: High-energy dissociative recombination rate co- 
efRcient calculated with the present MQDT treatment, com- 
pared with our results using calculated resonance curves and 
the O'Malley treatment (Ea. p2|l '). the theoretical results of 
Orel et a/. and the experimental results of Sundstrom et 
al. and Stromholm et al. [l^. Bottom: Contribution of 
each of the six sigma resonances to our total result using the 
O'Malley treatment. 



resonances give partial cross sections different from ours 
shown in Figure [SI 



V. DISCUSSION 

We have applied the theory of Ref. [l| to the calculation 
of DR rates for HeH+ + c^ in both the low- and high- 
energy regions. We have achieved very good agreement 
with experiment - at least as good as that of Orel et 
al. [l^ for the high-energy peak, and better than all prior 
calculations for the low-energy region. Thus, it is clear 
that the present methodology gives accurate results for 
both the indirect and direct mechanisms. 

Our treatment of the high-energy peak was perhaps the 
simplest possible: we used a quantum defect matrix con- 
stant with respect to both incident electron energy and 
internuclear radius. Such a simple treatment is applica- 
ble to much larger systems, as it leads to a sparse system 
of linear equations in the channel closing step, which is 
the rate limiting step for polyatomic DR calculations. 
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However, for larger systems it may be necessary to in- 
clude the energy dependence of the quantum defect in 
order to reproduce the relvant physics of DR or rovibra- 
tional autoionization. For instance, larger polyatomic 
systems such as NOj will in general support a greater 
number of shape resonances than smaller diatomic sys- 
tems, and these will impart energy-dependence to the 
quantum defect functions that is difficult to analytially 
remove. 

Methods suited to a MQDT treatment including the 
energy dependence of the fixed-nuclei quantum defect 
include those of Refs.[5l|, [s^. These calculations re- 
quire not the S-matrix but other entities such as the 
square root of the S-matrix or the quantum defect ma- 
trix. Therefore, such a calculation is more difficult owing 
to the necessity of following the branches of the quan- 
tum defect matrix across both R and E, for a large (> 1 
hartree) energy range. 

Even for the relatively small HeH+ system, the high- 



energy DR peak spans a large energy range and involves a 
large transfer of energy from the electronic to the nuclear 
degrees of freedom. Given these qualities, it is not a pri- 
ori obvious that our simple treatment would accurately 
reproduce the physical results. We pursued an energy- 
dependent treatment along the lines of [s^j , and this was 
yielding results similar to those presented in Fig|Sl but 
was not numerically stable, and would require improve- 
ment to yield publishable results. 
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